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Abstract—In this paper we present an efficient 
numerical technique, called the Characteristic Basis 
Function Method (CBFM), for the electromagnetic 
analysis of planar microstrip structures as well as for a 
class of three-dimensional scattering problems. In this 
method, the original problem geometry is segmented into 
smaller regions called blocks, and high-level basis 
functions are generated to represent the electromagnetic 
characteristics of these sections. These basis functions are 
referred to as the Characteristic Basis Functions (CBFs), 
and their use leads to a reduced matrix equation system. 
Since the method only requires the solution of relatively 
small-size matrix equation, both in the process of 
generating the CBFs and in the solution of the reduced 
matrix, the overall computational burden is reduced in 
terms of the memory requirement relative to the 
conventional Method of Moments (MoM) approach, and 
an acceleration of the solve time is also achieved. The 
Characteristic Basis Function Method is next combined 
with the Equivalent Medium Approach (EMA) for fast 
and efficient design of microstrip circuits etched on 
layered media. In particular, the developed EMA method 
substitutes the stratified environment with an equivalent 
“homogeneous” medium whose Dyadic Green's 
Functions (DGF's) can be evaluated analytically. The 
above technique yields reliable results and reduces the 
computational time in comparison with the conventional 
Method of Moments. In the second part of the paper, the 
CBFM will be applied to the solution of electromagnetic 
scattering problems involving plasmonic nano-rod array 
antennas operating in the Terahertz regime. In 
particular, the formulation necessary for the analysis of a 
single nano-rod is introduced and subsequently 
generalized to the case of a double-periodic array and to 
the case of a finite randomly tilted array. Several 
numerical examples which demonstrate the accuracy and 
the efficiency of the described procedures are included, 
both for the case of microstrip structures as well as for 
the nano-rod antennas. 


Index Terms—Method of Moments (MoM), 
Characteristic Basis Function Method (CBFM), Dyadic 
Green's Function (DGF's). 


I. INTRODUCTION 


Competition in the communication market has fueled the 
need to develop more and more efficient microwave 
Computer Aided Design (CAD) tools. One approach to 
accelerating the design process is to use a circuit simulator, 
as opposed to a full-wave field solver. In this approach, the 
microstrip discontinuities are modeled by lumped elements, 
and the system is analyzed by using a network-theory-based 
algorithm. While a network-theory type of circuit simulator 
requires a much smaller Central Processing Unit (CPU) time 
and memory than its full-wave electromagnetic counterpart, 
its accuracy degrades rapidly with an increase in the 
operating frequency, owing to the neglect of the parasitic 
coupling effects. Thus, it is almost always necessary to 
validate a design based on the network-theory-based circuit 
simulator, by using a full-wave electromagnetic solver, 
before the package can be sent out for manufacturing, to 
ensure that the device is "correct by design". 


One of the most widely employed full-wave algorithms is 
the Method of Moments (MoM) [1], which numerically 
solves the integral equations formulated by using the Dyadic 
Green’s Functions (DGF's) for the induced currents. Though 
numerically rigorous, and hence accurate, the MoM is 
computationally intensive, often requiring CPU time and 
memory that are orders of magnitude greater than those of 
conventional network-theory-based solvers. This is because 
the MoM matrix grows rapidly in size when the electrical 
dimensions of the problem become large, or when a fine 
mesh is used to model complex structures in order to 
guarantee a reasonable solution accuracy. In particular, if M 
is the total number of unknowns, the memory allocation and 
the CPU time increase as O(N?) and O(N’), respectively. 
This makes it difficult to incorporate the MoM solver into 
the design process, which often requires repeated 
modification and analysis of a circuit to achieve the desired 
performance. For the reasons indicated above, the use of full- 
wave electromagnetic simulations is often limited either to 
final design verification, or to a model development of the 
discontinuity for future use. 


In the past, many researchers have attempted to reduce the 
complexity and accelerate the computational speed of the 
MoM. For instance, the Fast-Multipole Method (FMM) [2] 
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and the Impedance Matrix Localization (IML) method [3] 
have been developed to expedite matrix-vector 
multiplications when using an iterative solver. Even though 
the FMM realizes a saving in the memory allocation 
procedure by storing only the near field interaction of the 
large matrix, and by employing an efficient spherical 
harmonic expansion technique for speeding up the solution 
of the linear equations, it is bounded by a discretization size 
ranging from 4/10 to 2/20, where / is the wavelength at the 
operating frequency. One consequence of this is that the 
MoM matrix grows rapidly as the geometry becomes 
electrically large. The IML transforms the source/testing 
basis functions into ones resembling traveling or standing 
waves, but its application is limited to sufficiently smooth 
objects. The wavelet transform [4] has also been used to 
yield sparse matrices that can be solved rapidly, though it is 
not as robust as one might desire. More recently, a novel 
iterative scheme called Maxwell and Markov method 
(MNM) [5]-[6] was proposed for frequency sweeping and 
successfully applied to a wide variety of problems. However, 
it is not as general as we desire. 


Another emerging approach for efficient MoM analysis of 
microstrip structures and scattering problems is based on the 
concept of segmentation or domain decomposition, and 
several techniques have been proposed to implement this 
concept. For instance, in [7], the modified diakoptic theory 
[8], originally proposed for antenna problems, has also been 
applied to microstrip structures, though its use has been 
relatively limited. The same is true for the diakoptic-theory- 
based Multilevel Moments Method (MMM) [9], which 
carries out an iterative basis-function refinement to solve 
passive planar structure problems. The Sub-domain 
Multilevel Approach (SMA), which utilizes the so-called 
Macro Basis Functions (MBFs) [10], is a novel technique for 
reducing the matrix size associated with large planar antenna 
array problems. A similar technique, called SFX (as before) 
[11]-[12], has been proposed to analyze microstrip structures 
and other types of antenna geometries in a numerically 
efficient manner by using Synthetic Functions (SFs). 


In this paper, we present the Characteristic Basis 
Function Method (CBFM) for efficient analysis of 
Monolithic Microwave Integrate Circuits (MMICs) Error! 
Reference source not found.-[20], as well as for a class of 
scattering problems involving three-dimensional objects 
[21]-[23]. In this method, the problem geometry is first 
segmented into sections, for which high-level basis 
functions, called the Characteristic Basis Functions (CBFs) 
are generated. These CBFs not only represent the unique 
electromagnetic characteristics of each section (Primary 
CBFs), but also include the interactions between different 
sections (Secondary CBFs), and account for parasitic 
coupling effects in a systematic and efficient manner. The 
CBFM leads to a reduced matrix, which is much smaller than 
the original one, and this obviates the need for iterative 
solution of problems requiring a large number of unknowns 
in their original MoM formulation. 

In the first part of the paper we discuss a CBFM-based 


technique which has been recently introduced for efficient 
simulation of microstrip circuits and printed antennas etched 


on layered media. The algorithm, referred to herein as the 
CBFM Equivalent Medium Approach (CBFM-EMA) [18]- 
[20], is well-suited for the fast and accurate calculation of the 
S-parameters of microstrip circuits. The CBFM-EMA 
replaces the stratified medium with a semi-infinite 
homogeneous grounded slab whose DGF's can be derived 
analytically. The proposed approach leads to a considerable 
time-saving when compared to the conventional approach, 
thanks to the analytical nature of the equivalent Green’s 
function, which is only comprised of a combination of the 
direct and reflected rays, the latter originating from the 
ground plane image. 


In the second part of the paper we show how the CBFM 
approach can be generalized with relative ease to the analysis 
of a class of three dimensional scattering problems. In 
particular, the adopted technique bypasses the computation 
of the Secondary CBFs, allowing us to generate the reduced 
matrix only by using Primary CBFs generated illuminating 
the block with a spectrum of plane waves. 


This paper is organized as follows. It begins by 
describing the main steps involved in the CBFM, and then 
presenting some numerical results which prove its accuracy 
and numerical efficiency. Section II presents a new efficient 
CBFM-based technique for the analysis of printed 
microwave circuits and antennas etched on layered media as 
well as the approach adopted in order to parallelize the 
CBFM algorithm. Section III extends the Characteristic 
Basis Function Method to the solution of electromagnetic 
scattering problems with particular emphasis on the analysis 
of periodic structures composed of plasmonic nano-rods. 
Finally some concluding remarks are summarized in Section 
IV. 


II. THE CHARACTERISTIC BASIS FUNCTION METHOD 


A. Formulation 


As mentioned in the previous section, the first step in the 
CBFM is to segment the original geometry into smaller 
regions, called sections, whose size are fractions of the 
original geometry, as for instance in the example shown in 
Fig. 1. 
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Fig. 1. Division of the overall region domain into N blocks. A couple of 
generic blocks involved in the Primary and Secondary CBF evaluation are 
colored in dark gray. 

Typically, two types of CBFs are defined for each 
section, namely, the Primary and the Secondary CBFs, 
though higher-order (for example, Tertiary) basis functions 
may also be included in the process, if necessary. The 
Primary CBFs are solutions used to represent the induced 
currents in the isolated sections, whereas the Secondary 
CBFs account for the field coupling between the sections. 
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The CBFs are generated by solving relatively small-size 
matrix equations, as compared to the original matrix in the 
MoM formulation. These matrices are diagonal blocks; 
hence, the time needed to generate these CBFs is usually 
reasonably small. In general, the level of field coupling 
between the sections is governed by such parameters as the 
geometry and the operating frequency, and we can take 
advantage of this fact by discarding some of the Secondary 
CBFs, in a dynamic manner, by using a thresholding scheme 
as reported in [14]. 


Next, the set of Primary and Secondary CBFs are 
employed as high-level basis and testing functions to 
generate a reduced matrix via the use of the Galerkin 
method. The resulting matrix is typically quite small and thus 
can be solved directly. 


The generation of the Primary CBFs starts by defining 
appropriate interfaces (inner ports) between adjacent blocks 
and considering each section in isolation. The number of 
Primary CBFs for each block is related to the number of 
interfaces associated with the section under consideration. In 
order to avoid an edge effect in the induced current 
distribution, the fictitious sources are shifted away slightly 
from the block edges by introducing a small extension. The 
matrix equation relating the voltage gap sources, applied at 
the block interfaces, and the Primary CBFs can be expressed 
as: 


Z y. 


.P. = 
=I]1 —,n —z,n n = 1,2, N; (1), 


where Z, is the impedance matrix corresponding to the i” 


h 


extended block; P,, is the n” current distribution 


associated with the 7” section; V, „is the n” source of block 


i; M is the block number; and N; is the number of Primary 
CBFs defined over the i” block. After solving the isolated 
problem, only the partitions of the CBFs entries belonging to 
the original blocks—and not to the extended ones—are 
retained. 


As mentioned earlier, the Secondary CBFs are needed to 
describe the electromagnetic coupling between the blocks. 
Specifically, a Secondary CBF represents the current 
distribution induced on a particular block due to the Primary 
CBFs residing on a different block. The Secondary CBFs are 
generated by solving the linear system of equations: 


i, j=1,2,,M i# J 

ZN, fines yÉ! 
Sp ee ee (2), 

m=1,---,N; 


where Z, is the block-diagonal matrix associated with the i” 
section, Z; is the off-diagonal matrix which accounts for the 


coupling effects between the block i and j, S,, is the n” 


in 
Secondary CBF associated with the i” section, N; is the 
number of Primary CBFs defined for the 7” block, and N;” is 


the total number of CBFs for the 7” block. 


The level of coupling between the various sections is 
determined by several parameters, such as the geometries of 
the blocks and the frequency of operation. We can therefore 
take a cue from the physics of the problem and discard the 
CBFs in a dynamic manner through a thresholding scheme. 
To implement this procedure, we apply the Singular Value 
Decomposition (SVD) algorithm and retain at most K 
linearly independent CBFs by choosing a threshold, which is 
typically on the order of 0.01-0.001. The aforementioned 
procedure enables us to reduce the number of employed 
Characteristic Basis Functions and, hence, the matrix size, 
while barely compromising the accuracy of the solution. 


The resulting current distribution J can then be 
expressed as: 


Ne Nee 

1 

J= ) Qi nd 4 oF; ) AM nJ 3 
Z at Ln+.n | M,n“M,n (3), 


where the œ; is the weight associated with the CBF, Jj i.e., 
the 7” CBF of the i” block after the truncation and SVD 
procedure. 


The final step in this process is the generation of the 
reduced matrix Z“ , whose dimension is (N x N"), where 


Nis the total number of CBFs, including both the Primaries 
and the Secondaries. The reduced matrix is generated by 
applying the Galerkin testing procedure, using the CBFs both 
as basis and testing functions (4): 

Z'a =V" 


(4), 
Z= ZJ, VSS = 


IS 


where Z is MoM impedance matrix, J is the matrix 


comprising all the selected CBFs and V is the vector 


associated with the analyzed input port. Typically, the 
dimension of the reduced matrix in the CBFM approach is 
much smaller than that needed in the conventional Moment 
Method formulation, and the reduced matrix can be dealt 
with directly via an Upper and Lower (LU) decomposition. 


B. Numerical results 


The first numerical example which illustrates the 
application of the CBFM is a band-pass filter, shown in 
Fig. 2. The current distribution induced on the filter is 
modeled by using Rao-Wilton-Glisson (RWGs) and rooftop 
basis functions. The total number of unknowns is 2390 and 
remains unchanged over the entire frequency range of 
interest, which spans from 5.0GHz to 11.0 GHz. The 
original problem is decomposed into two blocks, containing 
1198 and 1192 unknowns, respectively. The results for the 
current density distribution, computed along the middle of 
the input microstrip line of the first block at 8.0 GHz, and the 
filter S-parameters are displayed in Figs. 3-5, where the 
CBFM results are compared with those obtained by using the 
direct MoM solution (Direct Solution). 


We note that the proposed technique is able to accurately 
predict the frequency response of the filter even close to the 
device resonance frequencies. An SVD threshold of 0.01 is 
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employed and the relative number of selected CBFs is 
plotted in Fig. 6. 


We now estimate the accuracy of the algorithm by using 
the normalized (percent) error defined as: 


Rel. Error(%)=100- (5) 





where 1°? and 1,” represent the current density 


coefficients derived through the proposed technique and the 
conventional MoM, respectively. The relative error for the 
analyzed example is reported in Fig. 6. It is interesting to 
notice that the fluctuation in the relative error is determined 
by the total number of employed CBFS which is chosen in 
an arbitrary fashion. It is possible to lower the obtained 
relative error, which is already quite small, by increasing the 
number of CBFs. 
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Fig. 2. Analyzed band-pass filter. The filter is printed on a grounded lossless 
dielectric slab (¢. = 9.8, thickness ¢ = 15.0mil). 
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Fig. 3. Current distribution along the middle of the input line of the first 
block at 8.0 GHz. 

In order to provide a measure of the computational 
efficiency of the proposed approach, the Relative Time (RT) 
is defined as the ratio between the CPU time required by the 
CBFM in order to solve all the required linear systems of 
equations, and that needed by the conventional MoM. The 
RT for the analyzed example is 26.53%. 


For the second example we consider an 8x1 patch antenna 
array (see Fig. 7). The conventional approach to modeling 
this antenna requires 3416 unknowns, when a uniform 


rectangular grid and rooftop basis functions are employed. 
The problem is decomposed into 8 blocks, comprising of 427 
unknowns each. Since the array elements are isolated, the 
partitioning of the blocks is intuitively obvious. The 
principal-plane radiation patterns and the current distribution 
on the first patch antenna at 2.25 GHz are displayed in 
Figs. 8-11. 


It is evident that the proposed technique generates accurate 
results both for the radiation patterns as well as for the 
induced current distribution on the patch array antenna. The 
normalized percentage errors for the array problem is 
5.94x10°% when 64 high-level basis functions are 
employed. The RT for the array example is 1.93%. 
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Fig. 4. Return loss of the analyzed band-pass filter. 
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Fig. 5. Insertion loss of the analyzed band-pass filter. 


The last example we consider is that of a microstrip 
hairpin filter (see Fig. 12). For this geometry, the 
conventional MoM requires 2069 unknowns when using 
both RWG and rooftop basis functions. Furthermore, since 
the filter is composed of isolated elements, the segmentation 
of the geometry is intuitive and leads to a total of only two 
sections. 


The CBFM is implemented over the frequency band of 
interest, namely 11.5—13.0 GHz. The current distribution on 
the input microstrip line at 11.5 GHz, and the S-parameters 
derived from both the CBFM and the direct solution are 
compared in Figs. 13-15. The total number of CBFs and the 


4 
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correspondent relative error are plotted in Fig. 16. Note that 
very good agreement between the proposed approach and the 
direct solution is achieved over the entire frequency band 
analyzed. The RT for the hairpin filter example is 25.96%. 
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Fig. 6. Total number of CBFs and relative error for the band-pass filter 
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Fig. 7. Analyzed patch antenna array. The antenna is printed on a grounded 
lossless dielectric slab (& = 2.2, thickness ¢= 1.6mm). 
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Fig. 8. Magnitude of y-directed current distribution on the first patch 
antenna at 2.25 GHz obtained via the direct calculation. 


C. Parallelization of the CBFM 


One of the important attributes of the CBFM is that the 
algorithm is easily parallelizable. 


In this section we will outline the basic steps required for 
the parallelization and omit the details. Fig. 17 summarizes 
all the steps involved in the CBFM approach and their 
parallelization aspects are considered in this section. 


The basic concept in the CBFM is to decompose the 
complex structure under analysis into smaller subdomains 
for which a set of Primary and Secondary CBFs are 


generated. In order to maximize the parallel performance of 
the CBFM algorithm the same number of subdomains should 
be allocated to each processor for the purpose of load 
balancing. 
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Fig. 9. Magnitude of y-directed current distribution on the first patch 
antenna at 2.25 GHz obtained via the CBFM. 
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Fig. 10. H-plane radiation pattern obtained via the CBFM and the 
conventional MoM at 2.25 GHz. 
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Fig. 11. E-plane radiation pattern obtained via the CBFM and the 
conventional MoM at 2.25 GHz. 


If we consider an example where the total number of 
blocks is M, the MoM impedance matrix Z can be 
partitioned into sub-matrices as: 

5 
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Fig. 12. Analyzed hairpin filter. The filter is printed on a grounded lossless 
dielectric slab (s = 3.48, thickness t = 0.254mm). 
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Fig. 13. Current distribution along the middle of the input line of the first 
block at 11.5 GHz. 
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Fig. 14. Return loss of the hairpin filter. 


If the task of generating the CBFs for each subdomain is 
allocated to each process, the impedance matrix may be 
distributed among the processors as indicated in (7), where 
for the sake of illustration we have assumed that M=4, and 
we have two processes, identified by the numbers 0 and 1, 
respectively. 


Zi ] Z, 4 
=>? ~ Process ID =0 
Z _ E (7). 
= | S34 Process ID =1 
Ly 1 ° Z 4,4 


In (7) Z,, is the MoM sub-matrix describing the 


electromagnetic interaction between the blocks į and j. The 
primary excitation vector can be stored locally in the address 
space of each process. Distributing the impedance matrix 
among the processors in this manner balances the 
computational load from the point of view of the 
computational burden and memory usage. 
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Fig. 15. Insertion loss of the hairpin filter. 
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Fig. 16. Total number of CBFs and relative error for the hairpin filter 
example as a function of frequency. 

For the purpose of generating the Primary CBFs, no inter- 
processor communication is required, since each process 
contains all the required information to compute the CBFs. 
Generating the Secondary CBFs requires the knowledge of 
the Primary CBFs associated with all the domains. A 
collective MPI broadcasting command, MPI Allgatherv can 
be used to ensure that each process has all the required 
Primary CBFs in its address space. 
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Fig. 17. Overview of the CBFM approach. 


A low-rank reduced matrix equation can be obtained by 
using the Galerkin testing procedure with the generated 
CBFs as basis and testing functions. The described step is 
computationally expensive since it requires MK’ complex 
matrix-vector multiplications, where M is total number of 
blocks and K is an average number of CBFs per block. Each 
processor can however compute its own segment of the 
reduced matrix, significantly lowering the computational 
overhead associated with this step. For M=4, and for two 


rocesses. the reduced matrix Z? can be express as: 
9 s 














CBFM CBFM 
AS Zis | p ID=0 
TCBIM noo ne. Z CBFM FUERA 
£9) “94 
Ro — ~~ (8), 
= CBFM CBFM 
Zu Laa Process ID=1 
CBFM CBFM _ 
Lox Logs 
where: 
CBFM _ T 
ZG” =L ZL) ©), 


is the generic sub-block of the reduced impedance matrix, 
while J, is the matrix corresponding to the selected CBF for 


block i. Each segment of the reduced matrix can be gathered 
on a single process, which is designated as the "root" or the 
"master", employing the function MPI Gather. The 
remainder of the CBFM algorithm is performed sequentially 
on the "root" processor in order to avoid unnecessary 
overhead associated with the parallelization. 


In order to illustrate the numerical efficiency of the 
parallel CBFM algorithm, we analyze the case of the 24- 
patch antenna array shown in Fig. 18. The structure is 
modeled at the highest frequency of interest with a total of 
4064 rooftop basis functions, 1016 for each of the four 
blocks in which the antenna geometry is partitioned. The 
results for the radiation patterns of the antenna on the two 
principal planes at 2.5 GHz are shown in Figs. 19-20. 


The code has been run on 4 different CPUs and the results 
have been compared with a serial MoM approach (Direct 
Solution). Figs. 19 and 20 demonstrate that there is no loss in 
accuracy encountered in the process of parallelization. 


We have also tested the numerical efficiency of the 
parallel algorithm by subdividing the antenna geometry into 
24 blocks and progressively increasing the number of CPUs 
employed to solve the problem. 
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Fig. 18. Geometry of the analyzed patch antenna array. 
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Fig. 19. H-plane radiation pattern obtained via the Parallel CBFM and the 
conventional MoM at 2.5 GHz. 

The obtained normalized CPU Time is displayed in 
Fig. 21, together with the efficiency of the parallelization 
algorithm with respect to the ideal case where the solution 
time scales as //n, being n the number of allocated 
processors. Fig. 21 shows that the parallelized CBFM 
algorithm scales quite well indeed. 


D. Higher Order CBFM 


In the previous section, the main concepts of the CBFM 
have been summarized. It has been pointed out that, the first 
step in the CBFM is to decompose the domain of the original 
problem into smaller sections defined as "blocks" for which 
high-level basis functions are generated in order to 
accurately represent the electromagnetic characteristics of 
the analyzed device. To generate accurate results by using 
the conventional CBFM it is important to ensure that the 
inner ports be defined at locations where the effect of the 
higher-order modes is negligible such that the current 
distribution in the vicinity of these ports is expressible as a 
linear combination of a pair of forward and reflected waves 
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(see Fig. 22). To remove the constrains of this requirement 
on the domain decomposition process, which involves user 
intervention in defining the partitions, we have developed a 
new algorithm for efficient generation of the Primary CBFs. 
This new approach is particularly suited for the analysis of 
printed microwave circuits and antennas, and it enables us to 
generalize the decomposition procedure in conjunction with 
the use of Tertiary CBFs. 
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Fig. 20. E-plane radiation pattern obtained via the Parallel CBFM and the 
conventional MoM at 2.5 GHz. 

Furthermore, with the proposed approach, we can balance 
the computational load assigned to each block, maximizing 
the time performance of the parallel version of the CBFM 
algorithm in the process [16]. 
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Fig. 21. Normalized CPU time and parallelization efficiency of the CBFM 
code. 

Without loss of generality, let us consider the stop-band 
filter reported in Fig. 22. The filter is decomposed into two 
blocks, by partitioning the structure where the high-order 
modes are obviously present. In order to increase the number 
of linearly independent CBFs necessary to represent the 
induced current distribution in the vicinity of microstrip 
discontinuities, we employ the following voltage source 
expression for the internal ports (see Fig. 22): 


V =e" pl- N -=i 


n 


(10), 


where Aw is the distance between the center of the internal 
port and an arbitrary point of it (see Fig. 22, Dec. A). 
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Fig. 22. Geometry of the analyzed stop-band filter. The circuit is printed on 
a lossless dielectric slab with a relative dielectric constant s of 2.2 and 
thickness =0.794mm. 


The chosen expression for the internal-port voltage source 
is well suited for the analysis of printed microwave circuits 
and antennas etched on layered media. In particular, the 
factor N is intentionally overestimated (we have chosen 
N=15 for all the examples), though only the linearly 
independent CBFs are retained by employing the SVD 
algorithm. 


Furthermore, the proposed approach employs high-order 
CBFs, i.e., the Tertiary CBFs, which represent the current 
distribution induced on a particular block due to the 
Secondary CBFs residing on other blocks. We can then 
write: 


i, j=1,2,.,M iF] 
T 
7 — i n=l,- N; 
2 Lin a S jm r (11) 
m=], N} 


where /’,,, is the n” Tertiary CBF associated with the i” 


section N, N, and N; is the number of Primary CBFs, 
Secondary CBFs and Tertiary CBFs respectively, defined on 
the i” section. 


In all the analyzed cases, the use of Tertiary CBFs was 
found to be sufficient to obtain very accurate results; hence 
higher order basis functions beyond the Tertiary was not 
deemed necessary and were not considered. 


To demonstrate the accuracy and numerical efficiency of 
the proposed approach we begin by analyzing the microstrip 
filter problem shown in Fig. 22. The current distribution over 
the filter is modeled by using rooftop basis functions defined 
over a regular grid. The total number of unknowns is 636 and 
remains unchanged over the analyzed frequency range, 
which spans from 4.0 to 12.0 GHz. The original problem is 
partitioned into two blocks, with 328 and 308 unknowns, 
respectively. It is worthwhile pointing out that the filter is 
decomposed in close proximity of a microstrip discontinuity, 
where the higher-order modes obviously exist. 
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The magnitude and phase of the Sz; are presented in 
Fig. 23, where we compare the generalized CBFM algorithm 
proposed herein (Proposed Approach) with the conventional 
CBFM, as well as with some measurements, or with the 
Method of Moments when measured data are not available. 
Specifically, the filter has been analyzed with the 
conventional CBFM partitioning the structure where higher 
order modes (i) are present (CBFM Dec. A) and (ii) are not 
present (CBFM Dec. B). We note that the proposed 
technique is able to accurately predict the frequency response 
of the filter, even close to its resonant frequencies, while the 
conventional CBFM provides accurate results only when the 
circuit is partitioned where the effect of the higher-order 
modes is negligible. 


For the second example we analyze the 4x1 patch antenna 
array shown in Fig. 24. The conventional approach to 
modeling the array requires 1872 unknowns when rooftop 
basis functions are employed. 


The problem is decomposed into 2 blocks, comprising of 
1064 and 808 low-level basis functions, respectively. The 
magnitude and phase of the S;; and the S4; as well as the 
antenna radiation pattern at 7.75 GHz, are displayed in 
Figs. 25-28. 
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Fig. 23. S2; of the analyzed filter. 
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Fig. 24. Geometry of the analyzed antenna array. The circuit is printed on a 
lossless dielectric slab with a relative dielectric constant ¢=2.2 and 
thickness t = 0.794 mm. 


It is evident that the proposed technique generates results 
which compare well with the available measured 
data/conventional MoM. In particular, the developed 
algorithm yields accurate results even for the weak coupling 
coefficient S4; and the cross-polar component of the antenna 
pattern. 
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Fig. 25. Magnitude of the S;; and S4; for the antenna array. 


In order to further demonstrate the wider flexibility of the 
proposed approach, the results have been compared with the 
conventional CBFM by employing two decomposition 
strategies: arbitrary (Dec. A) or conventional (Dec. B). It is 
then evident that the partitioning strategy cannot be arbitrary 
if we employ the conventional CBFM. 


The accuracy of the algorithm is estimated by using the 
normalized (percentage) error defined in (5). 
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Fig. 26. Phase of the S;; and S4; for the antenna array. 


The total number of CBFs and the relative error for the 
analyzed examples are displayed in Fig. 29. It is worthwhile 
to point out that the relative error associated with the 
proposed scheme is below 1.0% over the entire frequency 
band, showing that the proposed scheme does not 
compromise its accuracy in order to gain robustness and 
flexibility. The fluctuations in the relative error plot can be 
attributed to the fact that the total number of CBFs selected 
by the SVD thresholding scheme varies as a function of the 
operating frequency. 
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In order to provide a measure of the computational 
efficiency of the proposed algorithm, the Relative Time (RT) 
criterion is employed. 


The RTs for the analyzed examples are 1.071 and 1.152, 
respectively. Hence, the time penalty (overhead with respect 
to the conventional CBFM) associated with the proposed 
algorithm is negligibly small, even though it gains us 
considerable freedom and flexibility in the implementation 
of the domain decomposition scheme, in comparison to those 
offered by the conventional CBFM, which itself has been 
proven to be faster than the conventional MoM. 
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Fig. 27. Antenna array radiation pattern on the ¢= 0° plane at 7.75 GHz. 


Patten(dB) 


Proposed Approach 
CBFM Dec. A 
CBFM Dec. B 
MoM 





“90 ~60 -30 0 30 60 90 
Theta(degs) 


Fig. 28. Antenna array radiation pattern on the ø= 90° plane at 7.75 GHz. 


MI. THE CBFM EQUIVALENT MEDIUM APPROACH 


A. CBFM-EMA formulation 


Rapid prototyping plays an extremely important role in 
the design of antennas and related planar circuits for wireless 
communications. Consequently, it becomes necessary to 
resort to the use of more sophisticated simulation techniques 
that are numerically rigorous, albeit computer-intensive. 

Without loss of generality, we consider the case where the 
background medium is stratified along the z-axis and is 
terminated by a perfectly conducting ground plane at the 
bottom (see Fig. 30). In the context of MoM, a circuit or an 
antenna printed on this type of stratified medium is typically 
analyzed by employing the Dyadic Green’s Functions whose 
spatial domain counterpart can be expressed in terms of 
Sommerfeld Integrals [24]-[33]. 


One of these approaches, namely the Mixed-Potential 
Integral Equation approach (MPIE) [33], utilizes both the 
electric and magnetic Dyadic Green’s Functions to express 
the electric and magnetic fields generated by the induced 
currents. 
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Fig. 29. Relative error and total number of CBFs for the analyzed examples. 


Then, for source and observation points located at z’ and z, 
respectively, the spatial-domain Green’s function f(p;z|z’) 
can be expressed in terms of the spectral counterpart 
F(kp;z\|z’) by using the following relationship: 








(0.0) 


2)= | rlp:z 


0 


= ie 0 (kek pak 





flo:z (12), 





where Jo is the Bessel function of the first kind of order 0, p 
is the horizontal distance in a cylindrical coordinate system 
and kp is the wavenumber in the o—direction. 


Microstrip 


wF 





Ground 


Ground 


Fig. 30. Grounded dielectric structure (left) and equivalent grounded 
homogeneous medium (right). 

The procedure described above, though general, and 
applicable to a multilayered structure with an arbitrary 
number of layers, is computationally intensive since the 
evaluation of Sommerfeld integrals is time-consuming [34]- 
[37]. 

Given this background, it is not surprising that a number 
of different approaches to speed-up the computation of SIs 
have been proposed. Among these, the Discrete Complex 
Image Method [25]-[28] is particularly appealing since it 
enables us to express the Sommerfeld Integrals analytically. 

By resorting to the Generalized Pencil of Function (GPOF) 
method [38] it is then possible to cast the spatial domain 
Green’s functions in closed form. Using the two-level 
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approach reported in [28], the spatial domain Green’s 


function can be expressed as: 
) > oS KiNn N, o Jki"on 
Pa an ———+ ) an e (1 
een T nel 2" = (13), 


2 2 2 2 
Nn=yP Oin Mn =P ay, 


and @1,2n/Q1,2, Ni and N: are the coefficients/exponents and 
the number of exponentials, respectively, obtained by using 
the GPOF method. 


Though the DCIM definitely represents an improvement 
over the purely numerical integration approach for 
evaluating the SI, the level of improvement can be 
inadequate, especially when the Green’s function has to be 
evaluated either at many different frequencies and/or for 
many different combinations of the vertical coordinates z and 
z’, Le., the observation and source locations. Furthermore, 
the DCIM often returns a large number of complex images 
for certain geometries, and results in a severe degradation of 
the numerical performance when this occurs. Although a 
number of techniques [39]-[44] have been recently proposed 
for the analysis of printed structures at a fixed frequency and 
planar microstrip devices at multiple frequencies, the 
difficulties faced in the evaluation of the SI have not yet 
been eliminated altogether. 


flo:z 





where: 


(14), 


In this section we propose an alternative approach to the 
DCIM referred to herein as the CBFM-Equivalent Medium 
Approach (CBFM-EMA) [18]. The underlying concept upon 
which the CBFM-EMA algorithm is based is relatively 
simple: construct the CBFs as well as the associated matrix 
equation for the unknown expansion coefficients for the 
semi-infinite and homogeneous grounded medium as shown 
in Fig. 30. 


It is worthwhile to point out that the Green’s function for 
the equivalent medium can be determined analytically, since 
it is simply comprised of the contributions of the direct and 
reflected rays from the source in the presence of the ground 
plane. Not entirely unexpectedly, the key step in this 
procedure is an accurate determination of the effective 
dielectric constant, which is a function of the ratio between 
the metallic trace width and the dielectric thicknesses. For a 
single-layer problem, the effective dielectric constant can be 
easily calculated analytically, as shown in [45]. However, the 
effective dielectric constant can be determined numerically 
by first solving the problem of an open-ended microstrip line 
residing on the multilayered substrate. The current 
distribution in the center region of the line, where it is free 
from edge-truncation effects, is expressed in terms of 
complex exponentials, via the use of the GPOF [38] 
technique, in order to extract the & value from the 
imaginary part of the propagation constant associated with 
the fundamental mode along the line. In the event the block 
under consideration is comprised of traces of different 
widths, the effective dielectric constant for the block is 
evaluated by using an averaging process. As expected, the 
proposed technique leads to a considerable time-saving as 


compared to the conventional approach, thanks to the 
analytical nature of the associated Green’s function. 


What is not so obvious though, and what has been 
discovered only through extensive numerical 
experimentations with a variety of different circuit and 
antenna configurations, is that the CBFM-EMA technique 
yields accurate results even when there are large variations in 
the widths of the metal traces of the considered problem 
geometry. Our experiments have demonstrated (see the 
numerical result section) that the CBFM-EMA results are 
found to be far superior, in terms of accuracy, when 
compared to the simple circuit-based approximations such as 
those generally implemented in commercial codes for rapid 
prototyping. Furthermore, typically, the differences between 
the results from the conventional MoM, and those derived by 
using either the Finite Element Method (FEM) or the Finite 
Difference Time Domain (FDTD) algorithms, are larger than 
the difference between the conventional MoM and the 
CBFM-EMA results for both the return loss or radiation 
characteristics. This leads us to observe that the CBFM- 
EMA is well suited for rapid prototyping of antennas and 
circuit configurations, because it offers considerable time- 
saving in comparison to conventional full-wave techniques, 
with little loss of accuracy introduced by the simplification. 


We now introduce yet another time-saving strategy, 
namely the Fast Matrix Generation (FMG) scheme, which is 
implemented in conjunction with the EMA [46]-[48]. Since 
the repeated numerical evaluation of the integrals involving 
the DGFs is the prime contributor to the matrix fill-time, the 
use of the FMG helps to speed up the evaluation of the MoM 
impedance matrix considerably. Thus, using this 
approximation, we can express the impedance matrix 
elements in a closed canonical form, and generate the MoM 
matrix in a fast and efficient manner while preserving its 
accuracy. 


Let us first consider the impedance matrix elements 
associated with a homogeneous problem in order to 
understand the physical interpretation of these elements. It is 
well-known that the impedance matrix relates the currents 
induced on the analyzed object to a suitable excitation 
source, as for instance a plane wave or a voltage source. The 
RWG or rooftop basis functions represent the flow of the 
current between the basis domains. For triangles or 
quadrilaterals of sufficiently small size, the radiated fields 
due to the currents flowing on a pair of subdomains may be 
approximated by the fields generated by a small dipole, once 
we go beyond a nominal distance (typically 0.154). Thus, 
using this approximation, we can express the impedance 
matrix elements in a closed canonical form, and generate the 
MoM matrix in a fast and efficient manner, with excellent 
accuracy. 


In order to associate an equivalent dipole moment to each 
of the low-level basis function, we need to identify the dipole 
moment location and its value, which is related to its 
magnitude and orientation. On the basis of extensive 
numerical experimentations we have found that the optimal 
location for the dipole moment center is the middle point of 
the common edge between the two triangles/quadrilaterals 
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that form the basis function (Fig. 31). The equivalent dipole 
moment, in the case of RWGs basis functions, can then be 
expressed as: 


M = {| fds = A +p) = (1 =r) (15), 
ror 


EF 


where p,“ and pp are the vectors between the free-vertex 
and the centroids of the subdomain 7," and Ty directed away 
and toward the vertex respectively while r,°° or r„® is the 
position vector of the centroids of T, or Tw. 


By using the previous dipole moment definition, the fields 
radiated by a low-level basis function can be expressed as 
follows: 





Test function 


T {e 
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Fig. 31. Electromagnetic interaction between two RWG basis functions via 

the FMG approach. 
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where /ol is the dipole moment value, k is the free-space 
wave number and 77 is the characteristic impedance of the 
surrounding medium. 


Furthermore, the previous FMG formulation can be easily 
generalized to the case of the equivalent semi-infinite 
medium (see Fig. 30), simply by using the image theory to 
account for the current distribution induced on the ground 
plane. The FMG procedure enables us to evaluate a 
considerable fraction of the impedance matrix entries 
analytically, which, in turn, reduces the fill-time of the MoM 
matrix significantly. In the next section, we will show how 
we combine the EMA algorithm, the FMG scheme and the 
Characteristic Basis Function Method for fast prototyping of 
microstrip circuits and antennas. 


B. CBFM-EMA numerical results 


To illustrate the CBFM-EMA procedure, we first analyze 
a four-stage, stepped-filter problem (see Fig. 32). The filter is 
modeled at the highest frequency by using 575 unknowns, 


and we leave this number unchanged over the entire 
frequency band of interest. 


The microstrip circuit is partitioned into 9 blocks, in a way 
such that the widths of the metal traces vary as little as 
possible within a block. The results for the S-parameters of 
the filter are shown in Figs. 33-34, where the CBFM-EMA 
results are compared with those obtained by using three 
independent approaches, namely the MoM, FEM and the 
Equivalent Circuit (EC) approaches. 


12.0 mm 
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Fig. 32. Geometry of the proposed four-stage stepped-filter. The filter is 
printed on a grounded lossless dielectric slab (s, = 4.0, thickness ¢= 1.0mm). 


Although the EC method is not as general as the others, it 
is still applicable to this type of geometry. It is apparent, 
however, that our technique generates results that are much 
more accurate than those obtained by using the EC 
approximation. 


We also note that the difference between the FEM and 
MoM results is larger than the corresponding difference 
between MoM and the proposed EMA, even near the 
resonances. 


This result is quite significant since the simulations with 
finite methods are usually much more time-consuming than 
they are with the conventional MoM. And yet, the EMA 
approach has been found to reduce the MoM-simulation time 
significantly, but yields results that are much closer to the 
MoM than does the corresponding FEM simulation of the 
same problem geometry. 


In Table I we list the values of the three resonance 
frequencies obtained by using the CBFM-EMA, the MoM, 
the FEM and the EC. 
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Fig. 33. S|; magnitude comparison for the stepped-filter problem. 


The relative differences between the CBFM-EMA, the 
FEM and the EC results from the MoM reference algorithm 
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obtained by using equation (17) are also included in the 
Table. 
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Fig. 34. S); magnitude comparison for the stepped-filter problem. 


Table I. Resonant frequency values for the four-stage stepped filter. 


fi h f Rel. fi Rel.f2 Rel. f 
EMA 5.26 5.55 5.83 0.38 0.36 0.17 
EC 5.03 5.4 5.77 4.01 2.35 1.2 
MoM 5.24 5.53 5.89 - - - 
FEM 5.22 5.64 5.94 0.38 1.99 1.71 


Next, in order to demonstrate the wider applicability of the 
proposed technique in comparison to conventional EC 
algorithms, we analyze three geometries that are difficult to 
handle by using the EC approach, though the CBFM-EMA 
can simulate them with relative ease. The EC approach will 
not be employed in the next examples. 


The first test example is a two-stub filter printed on a 
multilayered dielectric (see Fig. 35). The microstrip circuit is 
modeled by using 861 RWG basis functions, divided into 
two blocks, and analyzed with the CBFM-EMA. The 
corresponding S-parameters in the frequency range 2.0 - 
4.25 GHz are plotted in Figs. 36-37. We note that the 
proposed technique is able to accurately predict the stop- 
band of the filter, which ranges from 2.85 GHz to 3.70 GHz. 


For the next example, we address the problem of a patch 
antenna (see Fig. 38). The antenna is designed for operation 
at 2.40 GHz, and is printed on a substrate whose permittivity 
and thickness are 2.2 and 1.57 mm, respectively. A quarter- 
wavelength transformer is used to match the antenna input 
impedance to a 50 Q system. The total number of unknowns 
for the problem is 1081, when the mesh employed is 
uniform. The antenna is divided into 3 blocks and the 
corresponding S-parameters are calculated from 2.0 GHz to 
3.0 GHz (Fig. 39). It is worthwhile to note here that, 
although the CBFM-EMA is an approximate method and the 
analyzed geometry is characterized by large width variations 


of the metal traces, the results generated by the conventional 
MoM and by the proposed approach are still in very good 
agreement. 


In fact, the percentage difference in the predicted resonant 
frequency, computed by using (17), is only 0.83% for the 
CBFM-EMA and 3.35% for the FEM. 


The last problem analyzed is a patch antenna array shown 
in Fig. 40. The antenna is etched on a grounded, lossless, 
dielectric slab with & = 2.2 and thickness t = 1.57mm. The 
problem geometry is divided into 33 blocks and discretized 
by using a total of 6677 RWG basis functions. 


30.0 mm 





3.0 mm 
Fig. 35. Geometry of the two-stub microstrip problem. The filter is printed 
on a grounded lossless double layer dielectric slab (from the bottom 
&) = 2.1, &2 = 12.5, thickness ¢; = 0.7mm, t; = 0.3mm). 
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Fig. 36. Sıı magnitude comparison for the two-stub filter problem. 


The principal-plane radiation patterns and the current 
distribution along the input feed-line at 2.5GHz are 
displayed in Figs. 41-43. It is evident that the proposed 
technique generates accurate results both for the radiation 
patterns and for the current distribution of the patch array 
antenna. 


In order to provide a measure of the computational 
efficiency of the proposed approach, we define the Relative 
Time as the ratio between the CPU time required by the 
CBFM-EMA and that needed by the conventional MoM. The 
RTs for the four examples investigated are 20.75%, 37.9%, 
37.37% and 30.25%, respectively. In comparison, the RT 
between FEM and the conventional MoM is higher, and the 
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accuracy of the Finite Method is usually lower than that of 
MoM for this type of problems. 
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Fig. 37. So; magnitude comparison for the two-stub filter problem. 
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Fig. 38. Geometry of the analyzed patch antenna. 


IV. CBFM ANALYSIS OF NANO-ROD ANTENNA ARRAYS 


A. Introduction 


In this section, we will extend the previously described 
CBFM for the analysis of a class of scattering problems [21], 
[22], [49]-[52]. In particular, we will focus on the 
development of an efficient scheme to accurately analyze 
finite and doubly periodic arrays of plasmonic nano-rods 
illuminated by an obliquely incident plane wave. The 
problem is formulated by using integral equations and the 
Method of Moments in conjunction with the CBFM to 
significantly reduce the global number of unknowns. 


Arrays of plasmonic antennas demonstrate interesting 
physical attributes and find a wide range of applications in 
photonics such as light-trapping arrays of plasmonic nano- 
particles [53], arrays of gold nano-rods loaded with poly (3- 
hexylthiophene) (P3HT) [54], [55], optical waveguides 
realized with one-dimensional arrays of plasmonic nano- 
particles [56]. Furthermore, arrays of rectangular gold nano- 
rod antennas surrounded by an insulating SiO2 region can be 
employed as Schottky diodes [57] and nano-optical Yagi- 
Uda antennas, comprising of plasmonic nano-rods that can 
both enhance and direct the interaction of a single quantum 
emitter with electromagnetic fields [58]. 
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Fig.39. Sıı magnitude comparison for the patch antenna problem. 
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Fig. 40. Geometry of the patch antenna array under analysis. 
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Fig. 41. H-plane radiation pattern obtained via the CBFM-EMA, the 
conventional MoM and the conventional FEM at 2.5 GHz. 

We face two main challenges when attempting to solve the 
problem of the array of plasmonic nano-rods: first, how to 
model the elements (especially thin elements with small 
features) in the unit-cell by using as few unknowns as 
possible; and, second, how to address the issue of slow 
convergence of the Periodic Green’s Function (PGF). To 
address the first challenge, we employ Macro Basis- 
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Functions (MBFs) [50], [59], that are high-level basis 
functions, to model the induced current with only a few 
unknowns. 
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Fig. 42. E-plane radiation pattern obtained via the CBFM-EMA, the 

conventional MoM and the conventional FEM at 2.5 GHz. 
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Fig. 43. Comparison of the current distribution along the middle of the input 
feed-line at 2.5 GHz. 


Without loss of generality, let us consider a two- 
dimensional periodic array of plasmonic nano-rods lying in 
the x-y plane, and illuminated by an arbitrarily incident plane 
wave, given by: 


(18), 


where kinc 1s the free space wave vector and r is the position 
vector, as shown in Fig. 44. 


[hie = Ege Ene 2 


The unit-cell is assumed to be rectangular with 
periodicities along the x- and y-directions equal to Ax and 4y, 
respectively. 


The isolated element in the unit-cell is a nano-rod 
composed by a plasmonic material, which can be 
characterized at optical frequencies by the Drude model as: 


E, = Em — fo NEC — if) (19), 





Fig. 44. Nano-rod illuminated by a plane wave with parallel polarization. 


where fp is the plasma frequency, and fa is the damping 
frequency [60]. 


In the analysis of periodic structures, the scattered field 
can be expressed a summation of Floquet harmonics [1], 
comprising of a combination of propagating (belonging to 
visible range) and evanescent modes (invisible range). We 
assume that the dimensions of the spatial lattice, i.e., Ax and 
Ay, are such that only a single Floquet mode exists in the 
visible region. In particular, we can characterize a periodic 
structure either by computing its dispersion diagram, which 
shows the range of frequencies within which there are no 
propagating modes, or by deriving the reflection coefficient 
characteristics over a range of frequencies. There are two 
advantages to presenting the results for the reflection 
coefficient, which is the focus of our work. First, one can 
provide the information on phase and polarization, and 
second, information on evanescent modes can also be made 
available. We will now describe the details of our modeling 
procedure in the sections below. 


B. Solving the isolated element problem 


We now proceed to briefly review the technique for 
analyzing a plasmonic nano-rod, which has been presented in 
detail in [50]. Fig. 44 shows the nano-rod, which is oriented 
along the y-axis and illuminated by a plane wave with 
parallel polarization. 


To formulate the above antenna scattering problem, we 
employ the volume equivalence theorem [61] to replace the 
rod with a volume polarization current J radiating the 
scattered field E°“: 


jole, S E ay ag : a _ (20), 


i 
A 
where A is the nano-rod cross section area. 


We note that in contrast to the case of the PEC rod which 
only requires longitudinal currents, the volume polarization 
current can have both longitudinal as well as transverse 
components. Since the nano-rod is electrically thin, we can 
assume that the total polarization current vector Z is 
uniformly distributed across the antenna cross-section. 
Furthermore, we assume that the polarization current J exists 
along the axis of the rod. 


Our next step is to express / in terms of N Macro Basis 
Functions as depicted in Fig. 45. The w-component of Z is 


written as: 
N-I 
ly = >, itt (cr) (21), 
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where t, is the n MBF with axis along a and centered in a: 


D (a) = sinlko (H -|a —Ay ) 


(22). 





Fig. 45. Discretization of the polarization current by piecewise sinusoidal 
Macro Basis Functions. 


By Galerkin testing on the surface of the nano-rod, one 
can derive an accurate, fast and singularity-free 
computational paradigm to solve the problem at hand and 
obtain the unknown current. The current induced on the axis 
of the nano-rod can be obtained by solving: 


(23), 


[MBF _ 7 MBF MBF __y MBF 
G =res 


where /”"" is a column matrix that includes the coefficients 


of each MBFs, and 


Ym tl 

Zar | = fin (an dy (24), 
Yml 
Ym tl 

yar |= [ER Olmos, 
-1 


and tm(y) is the m” piecewise sinusoidal MBF. As depicted in 
Fig. 45, the first and last MBFs are half, while the rest are 
full piecewise sinusoidal functions. 


The length of the domain of the half MBF is H, and ym=- 
L/2+(m-1)H is the location of the cusp of each MBF. The 
electric field radiated by the n” MBF is denoted by G,”2" in 
(24) and is defined as: 


(GPF -Right (y ~y, ) n=] 
Gr = Gan W- Yn) n=N mF. (26), 
GPA Cy- y) otherwise 


where GMP. RighLeñ is the electric field, radiated by the 
right/left half-sinusoidal basis, depicted in Fig. 46, which can 
be expressed in a closed form as indicated in (27). 


In particular, G”?"!" represents the summation of the 


electric field radiated by the left and right half- sinusoidal 
MBFs. The observation point, for which the impedance 
matrix is constructed, is located at (x, y, z) = (0, y, a), on the 
surface of the nano-rod. 


~jn Paes 
ars ee 


oTkoRe 


Ri Ro 


Gy MBF Right 





o- jkoRo 
-ysin(ko H)—— |, 
ko Rò 


(27), 


[eosar + jZ sinat) 
0 


where 
Ri = x? +(y-HY pr 


Ro =x? +y? Ja 


Finally, the impedance matrix element can be obtained 
from: 


(28). 


Sif m=n#(l,Niypr) 
ptp [m-n =| 





z -j1 
Z ee LE 
0 otherwise 
where s f and s „+ are expressed as: 
sin(2k oH) 
S ¢ = H| 1- —— 
if Te (30), 
s, -Zf sinko) _ (49H) 
Wh I| koH 0 (31). 





Fig. 46. Full triangular Macro Basis Function including right and left half- 
MBEs. 


We now proceed to analyze a plasmonic nano-rod antenna 
which is made of silver, which can be described as a Drude 
medium at optical wavelengths in terms of the following 
parameters: 


Ero =5, fp =2175 THz, fg =4.35 THz (32). 
We assume that its length is Z=300nm, and that it has a 
radius of a=10nm. 


The z-component of the scattered electric field is 
computed at p =500nm over a range of optical frequencies. 
The scattered field is ¢symmetric for the normally incident 
wave. It is found that numerical convergence is achieved 
with only five MBFs. The resonant frequency obtained by 
using the MBF formulation, at which the magnitude of 
scattered electric field is maximum, is found to be 
fres=187THz (1.6um), which is almost the same as predicted 
by the FDTD analysis (Fig. 47). 


The computational time required to derive the scattering 
characteristics over the entire frequency band by using the 
MBF method is ~5s for 500 frequency samples when the 
simulation is performed on an Intel Core(TM)2 Duo at 
2.2GHz CPU frequency with 4GB of RAM, while the 
corresponding time is on the order of 4h for the FDTD 
method. In fact, the generation of surface plasmon mode on 
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the curve boundary of the nano-rod requires the use a fine 
mesh in the FDTD to obtain accurate results. We see that a 
reduction of the computational time in the range of three 
orders of magnitude with respect to the FDTD method is 
achieved by the proposed approach. 


1 
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0.8) R | e FDTD 
| k 
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Fig. 47. Frequency response of the plasmonic nano-rod illuminated by a 
normal incident vertically polarized plane wave. The scattered field is 
computed by using 5 MBFs (solid blue line) and FDTD (circular red 
marker). 

We now turn to a second example, for which the angle of 
incidence is @ =76.7° and g--45°. We assume that the 
operating wavelength is 1.5um (4-200THz) which is close to 
the plasmonic resonance. In contrast to the previous 
example, we need to consider the transverse components of 
the polarization current, in addition to the longitudinal one. 
Five MBFs are found to be sufficient to achieve 
convergence, for each of the x-, y-, and z-components of the 
polarization current; hence, the size of the matrix equation is 
15x15. 


The scattered near-fields along the x-axis for 
0.05A<x<0.5A are shown in Fig. 48 and are compared with 
the corresponding FDTD-generated data. It is evident that 
very good agreement is achieved between the results derived 
by using the two methods. Thus, the MBF method achieves a 
considerable time-saving without compromising the 
accuracy. 


We now turn to the far field radiation patterns of the 
plasmonic rod in the x-z plane, shown in Fig. 49, that are 
computed at r=10/A. Very good agreement is achieved 
between the results from the proposed approach and the 
Finite Difference Time Domain Method, however, as we 
have mentioned earlier, the MBF approach is much faster, 
since it takes only 0.2s while the FDTD requires 4h to 
generate the same results. 

C. Doubly Periodic Array of Plasmonic Nano-rods 

In this section we will generalize the previously described 
approach to the analysis of doubly periodic nano-rod arrays. 
In particular, the coupling effects between the elements will 
be taken into account by using the CBFM. 

The CBFs for the nano-rod are generated by illuminating 
the block under analysis with a spectrum of plane waves 
[PW;, PW ...], as shown in Fig. 50, with both parallel and 
perpendicular polarizations. 

Let us denote as J,(r‘) the induced current distribution 
associated with the n” plane wave. Next, we construct the 
matrix J: 


eae ee (33). 


where N, is the number of employed plane waves. 
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Fig. 48. Magnitude of the E, component along the z-axis. 
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Fig. 49. Normalized radiation pattern (E, component) in x-z plane. 
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Fig. 50. An arbitrary object illuminated by N, incident plane waves with 
different polarizations to obtain the CBFs. 


Finally, the CBFs are generated by applying the Singular 
Value Decomposition factorization [18], as follows: 


J=U- 2V" (34), 


where U and V are the singular vector matrices while the 


main diagonal of x contains the singular values. Not all the 


determined singular vectors are significant. In practice only 
the vectors associated singular values above a certain 
threshold, after being normalized with respect to the biggest 
one, are retained. The selected singular vectors are then 
employed as CBFs. In all the reported results, the threshold 
has chosen to be equal to 0.001. The first two selected CBFs 
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at f = 260 THz for the nano-rod problem are reported in 
Fig. 51. 


We now proceed to solve a periodic array problem by 
utilizing the previously generated CBFs. In particular, we 
avoid to use the periodic Green’s function approach, which is 
often characterized by slow convergence issues, developing 
an alternate strategy which is considerably more efficient 
[62]. 





CBFs, Magnitude, f=260 THz 
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Fig. 51 Magnitude of CBFs for plasmonic nano-rod at /=260 THz. 


By defining as Ax and A, the spatial period along the x and 
y directions [1], the pg-plasmonic antenna located at 
(x,V)=(pAx, gAy) can be represented by an induced current 
Tng(v)=I epg. Where I(y) is the unknown induced current 
distribution along the axis of the center element and e,, is the 
inter-element phasing factor [63] introduced by the incident 
plane wave: 
epg = exp- iNeine pA, + ky qh y | (35). 
By imposing the polarization current condition (20) and by 
expressing /(y) in terms of the selected CBFs, the scattered 
electric field, £,°°”, can be written as: 


gseat _ Ne pCBF ae GCBEn 
y 5 n e pqV pq 
n=1 p,q =—2 


where G,,“?""" is the y-component of electric field radiated 
by the n” CBF on the (p,g)” element. 


Finally, the following matrix equation can be derived by 
applying the Galerkin testing procedure [61]: 


po _ zCBF CBF _ -y CBF 


(36), 


—TSres (37), 
where: 
L/2 
MBF | _ 7 CBF .n 
z4 be 0) Sao dy (38), 
-L/2 
l L/2 
MBF -J1 
É | E ee | tmin(v)dy 89), 
=res Jmn ma kole, -1) ae 
L/2 
CBF j 
v ba [E Omno (40). 
—-L/2 


It is worthwhile to notice that the impedance matrix 
element (38) contains a doubly-infinite series, which 
represents the contributions of the array elements. A possible 
approach to evaluating the result of the impedance matrix 


element is to use the concept of progressively expanding 
rings, as illustrated in Fig. 52. 
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Fig. 52. Ring definition employed in the context of convergence of the 
Galerkin integrals in the spatial domain. 


In particular, at stage M, we need to solve a finite array 
problem composed by (2M-1) x(2M-1) elements. It is evident 
that the problem will rapidly become time consuming when 
the number of rings is large. The reason for such slow 
convergence in the spatial domain can be ascribed to the 
spherical wave terms in field scattered by each nano-rod. 


In Figs. 53-54 the unknown coefficients associated with 
the first and third CBF are reported as a function of the ring 
number M. 


It is evident, however, that after a few stages, the maxima 
and minima of the CBF coefficients appear at the same 
location of the "ring" axis, implying that the shape of the 
final current distribution will not change and only its level 
will fluctuate by increasing the number of rings. We can take 
advantage of this fact by solving a relatively small-size 
matrix equation to find the shape of the CBF-combination. 


Now we will discuss a computationally efficient approach 
to evaluate the unknowns weight associated with the each 
CBF. The macro-CBF, /(y), can be written as 


LW) = CL shape O) (41), 


where C is the unknown weight coefficient. By applying the 
Galerkin testing procedure to the polarization current, the 
level of the macro-CBF can be expressed as: 


{ CBM) j CBF 
ar a ae (42). 
{ — ( CBF _ ee iil 
— — G 


=res 
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Fig. 53. Magnitude of the coefficients of the first CBF as function of the 

ring index for normal incidence. 
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Fig. 54. Magnitude of the coefficients of the third CBF as function of the 
ring index for normal incidence. 


We note that, the level of macro-CBF depends on the 
coefficients of CBFs after M rings. 


To address the problem of evaluating (42), we switch to 
the wave-number domain by recurring to the Parseval's 
theorem which provides: 


+00 +00 +00 
| | | effa Fy (ky) ES” dk dky dk, 
CBF —00 —00 —90 
[al F 
G m,n (27)° ( ), 


where the symbol tilde, ~, denotes the Fourier Transform 
(FT). The Fourier transform of the CBF current distribution 
appearing in (43) can be expressed as: 





2 nee 
~CBEn;,,_ 40° F p q 
TOEKA) D EEA K) (44 
á p.q =—2 
where: 
20 
kP = p4 pinc 
x OK, x 
i= 2m _ inc (a 
y-4 Ay y 


while the FT of the corresponding dyadic Green's function 
element is defined as: 


ee 
G 4) ina" 
= ko k-k? 





(46). 


By substituting the scattered field expression in (43), we 
obtain: 


Y nD REDS ED -e lind) 
_ P,]=— 9 
mn 7 27k AyAy, 





car 


G (47), 


where /nt(k,p) is reported in (48). 


In this paper we assume that the spatial periodicities of the 
lattice are smaller than half wavelength at the operating 
frequency implying that only the main Floquet mode is 
allowed to propagate. 


+00 Jkr 


Intlk ») = | a - 
EEF 

















—o0 
-fe-l 
= Je 2 2 
A s ky > kp 
P ky -kp (48). 
= me" P? x 
-| k3 -ko |IzI 
z oi 
ke -k2 
pP 0 





By considering one visible mode only, the generic 
impedance matrix element can be expressed as: 


, pec\ Nt 
cB | = tim yy ye cer f (Ni) 
G dnn G nn 


where N, is the truncation number for the infinite series 
embedded in the spectral representation of the matrix. 


(49), 


We now proceed to the analysis of a silver nano-rod array 
with length L=150nm and radius a=7.5nm. In order to prove 
the computational efficiency which characterizes the 
proposed approach, we compute the (1,1)-element of the 
reduced impedance matrix by using the spatial and spectral 
techniques. The results are shown in Figs. 55 and 56, for 
normal and oblique incidence, respectively, where on the 
horizontal axis the series truncation number is reported. It is 
evident that the convergence in the spatial domain is slow 
and not monotonic; while the spectral scheme converges fast 
and without oscillation. Both for the normal and oblique 
incidence case only 20 elements in the series are sufficient to 
guarantee convergence of the final results. The described 
problem has been run on a 2.2 GHz Intel Core(TM)2. The 
reflection and transmission coefficients for both normal and 
oblique incidence are plotted in Figs. 57-58 and compared 
with the corresponding results from the full wave FDTD 
simulation. To achieve an accurate result with the FDTD 
modeling, fine meshing must be used for the thin cross 
section of the rod, leading to a simulation time of about 10h 
for the normal incidence case. For oblique excitation, we use 
Sin/Cos single frequency calculation [64] in FDTD with 
coarse meshing. The calculation is performed at 24 points, 
requiring approximately 5 days to complete. On the contrary, 
the CBFM is considerably faster and it takes only about 30s 
for 73 frequency samples to generate accurate results. It is 
finally interesting to remark that the CBFM is a direct 
approach and can, differently from the FDTD method, 
efficiently handle multiple right hand side type of problems. 


D. Modeling large non-uniform optical antenna arrays for 
meta-surface applications 


Metamaterials that are made of sub-wavelength 
composites (“meta-atoms’’) have gathered a lot of attention 
in recent years, because of their unique electromagnetic 
properties that are unattainable in nature, with promise of 
extraordinary functionalities and applications, such as 
negative refraction and invisibilities. 
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Fig. 55. Magnitude of the (1,1)-element of the reduced matrix at 4-400 THz 
computed by the spatial and the spectral methods for normal incidence. 
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Fig. 56. Magnitude of the (1,1)-element of the reduced matrix at 400 THz 
computed by the spatial and the spectral methods for oblique incidence. 
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Fig. 57. Reflection and transmission coefficients (in percentage) for normal 
incidence. 
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One of the most interesting meta-structures is based on the 
emerging concepts of optical antennas and a special class of 
meta-materials are the 2D meta-materials, also called meta- 
surfaces. 


Simulation of nano-antennas for  meta-surfaces, 
comprising of dispersive metallic nano structures and 
dielectric environments, plays an important role in the design 
of meta-surface-based devices. 
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Fig. 58. Reflection and transmission coefficients (in percentage) for oblique 
incidence. 


A whole host of methods, such as the Finite Element 
Method, Finite Difference Time Domain algorithm, Method 
of Moments, and the Boundary Element Method (BEM), 
have been employed for such simulations, after they have 
been tailored for plasmonic applications. Even so, modeling 
large nano structures still poses great challenges, because 
unlike uniform phased arrays, or Frequency Selective 
Surfaces (FSSs), meta-surfaces comprise nano-antennas that 
could in general be different in size, shape, and orientation. 
For such non-periodic geometries, we can no longer solve 
the problem at hand by considering only a single unit cell. 
Conventional numerical methods become inefficient when 
analyzing large arrays and the size they can handle is often 
limited by the available CPU time and memory. 


The objective of this section is to apply the CBFM 
previously developed, to the analysis of finite, non-uniform, 
and large plasmonic nano-rod arrays, which show promise in 
designing ultra-thin optical devices [52]. We choose 
plasmonic nano-rods for the array element, with spatially 
varying angles of orientation, which can generate interfacial 
phase change for circularly polarized waves. 


When a circularly polarized wave illuminates a nano-rod 
whose axis has an angle y with respect to the x axis, the 
scattered field can be expressed, for small incident and 
observation angles as: 

po E(t ee js2y | (50), 
where s = +1 denotes the helicity of incident wave (“+” 
corresponds to right-handed (RCP) while “-” corresponds to 


left-handed CPs (LCP)) while ~ = Xx JY represents the 
direction of the polarization. 


We can achieve any desired phase distribution by adjusting 
the orientation of the nano-rods, and this property is useful 
for designing novel reflection and/or refraction-based optical 
devices. 


For the first example, we consider an optical vortex 
formed by plasmonic nano-rods. In order to realize a helical 
wave front, we introduce an azimuthal phase variation, i.e., 
exp(-jl@), on the meta-surface, where / is the topological 
quantum number and g-2y represents the azimuthal 
coordinate. Fig. 59 shows the analyzed optical vortex with /= 
-! and the corresponding in-plane phase retardation 
introduced by the meta-surface. The array consists of 40x40 
nano-rods with length of 275nm and radius of 10nm that are 
placed equally spaced. 

The separation distance between the centers of the 
neighboring elements is 500nm. A RCP Gaussian beam is 
normally incident from the top, whose waist is located on the 
meta-surface; its radius is 10nm, and its wavelength is 
A=1.55um. In order to analyze the array we use the CBFM 
with 40x40=1600 blocks. 


We observe the scattered field for the cross-polarized 
(LCP) components (see Fig. 60). The observation plane is 
120um*x120um large, which is parallel to the meta-surface 
and is located at a distance of 400um below the array. 
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Fig. 59. Schematic of an optical vortex and additional interfacial phase 
dependence introduced by the optical vortex. 
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Fig. 60. Far field intensity pattern of an optical vortex carrying topological 
charge /=-1, calculated by the CBFM and the conventional MoM. 


The calculated radiation pattern shows good accuracy 
between the CBFM and the conventional MoM results. For 
this example, the generation of the CBFs takes a total time of 
0.62s. The CPU time required to solve the resulting linear 
system of equations and the relative allocated memory are 
4.5s and 20.48MB, respectively, for the CBFM while these 
figures rise to 188.6s and 2.048GB when the conventional 
MoM is used. 


For the next example, we consider an array with a 
sunflower Fibonacci pattern, as shown in Fig. 61. This meta- 
surface was created by arranging the nano-rods such that the 
n” nano-rod is located at (On, G,), with the coordinates of the 
location governed by the Fermat’s spiral, as follows: 


p, =LvNn 


(51), 
of =ng 
where g=27(/-1/9); @ is the golden ratio; Lo=250nm is a 
scaling constant. The antenna is composed by 1600 nano- 
rods with their orientations tangent to the spiral. 


Moreover, the radius of the nano-rods is 10nm while their 
lengths vary randomly from 265nm to 285nm. The meta- 
surface, which introduces an azimuthal phase gradient 
creating an optical vortex with /=2, is illuminated by an LCP 
Gaussian beam, whose waist radius is wo=l0um. The 
wavelength for this simulation is A=1.55um. The problem is 
divided into 1600 blocks, with each block containing only a 
single nano-rod. 


The generation of the CBFs takes 133.4 s for 1600 blocks, 
when using a server with 12 processors. The far field of the 
opposite cross-polarized wave are computed over a 
200umx200um large plane, at a distance of 400um below 
the meta-surface (see Fig. 62). Good agreement is achieved 
between far field results derived by using the CBFM and the 
MoM. 
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Fig. 61. Schematic of a Fibonacci spiral nanorod array creating an optical 
vortex. 
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Fig. 62. Far field intensity pattern distribution of the Fibonacci optical 
vortex. 


For the third example, we consider an optical axicon 
realized with straight nano-rods which can be used to 
generate an approximation of a Bessel beam with a 
collimated Gaussian beam input (Fig. 63). In order to form a 
conical wave-front of light that travels through the meta- 
surface, the phase delay has to increase linearly with the 
distance from the center. For a given radius p on the meta- 
surface, the phase retardation must satisfy: 


2m, 
~ = — psinO (52), 
A 
where @ = tan'(R,/DOF), R, is the radius of the meta- 
surface, and DOF is the depth of focus of the axicon. 


We analyze an axicon with a radius of R,=20um, 
consisting of 5024 nano-rods with varying orientations. Each 
nano-rod is identical in length, which is 275nm, and has a 
radius of 10nm. The axicon is illuminated by a RCP 
Gaussian beam, whose waist radius is wo=40um; and its 
wavelength is A=1.55um. 


Figure 64 plots the computed field intensity of the 
designed axicon, which is located at z=0, showing that light 
does not spread out as it propagates within the axicon’s depth 
of focus. 


“TU 
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Fig. 63. Schematic showing the design of the axicon. 
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Also, the beam forms rings of increasing diameters over 
the radial distance, while preserving their thicknesses, thus 
closely replicating the properties of a Bessel beam. 


Finally, to demonstrate the ability of the proposed method 
to model large plasmonic arrays, we analyze a flat lens 
comprised 31428 elements, as depicted in Fig. 65. In order to 
form a focus at a distance such that all light arrives there in 
phase, additional phase retardation should be introduced to 
compensate for the phase lead introduced by the disparity in 
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Fig. 64. Calculations results by using the CBFM. The results showing the (a) 
xz (longitudinal) and (b) xy (lateral) cross sections of the far-field intensity 
distributions. (c) Plots of the intensity along the x axis on y = 0, z=-33um 
plane. 





Fig. 65. Schematic of the design of the lens with a focus of 200um. 


For a focal length f, the phase dependence ø as a function 
of the radius p can be expressed as: 


p= Nr - 4) 


The lens, whose diameter is 100um, is located at the z=0 
plane, and its designed focal length is 7-200 um. It is 
illuminated by an RCP plane wave whose wavelength is 
A=1.55um. Figure 66 shows the results calculated by using 
the CBFM. As we see from Fig. 66 light is confined to the 
focal region for the normal incidence case, as desired. 

It is crucial to remark that, once the impedance matrix has 
been filled, we can solve for multiple excitations by using 
the same matrix, without generating a new one for each 
incident angle, which, in turn, results in significant time 


(53). 


saving. Figure 66(b) plots the field intensity of the xz cut- 
plane of the meta-lens at y=0, for the oblique incident angle 
0”"°=5°. We observe that the focal beam tilts as desired. 
Figs. 66(c) and (d) show the corresponding xy cross-section 
of the field distribution at the spots. 

The conventional MoM requires 395GB of memory while 
only 3.95GB are needed by the CBFM. Furthermore, the 
solve time is on the order of 6 mins on a single processor, 
which is much faster than that required by the MoM. 
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Fig. 66. Calculation results for the lens by using the CBFM. Longitudinal 
cross section of the field distribution for the (a) normal incidence and (b) 


oblique incidence, respectively. (c) and (d) transverse cross sections of the 
field distribution at the spots corresponding to (a) and (b), respectively. 
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V. CONCLUSIONS 


In this paper, we have presented a numerically efficient 
algorithm—the Characteristic Basis Function Method—to 
analyze planar microstrip structures by combining the 
concepts of domain decomposition and high-level basis 
functions. The generation of the characteristic basis functions 
involves the LU decomposition of relatively small-size 
matrices, whose dimensions are a fraction of those of the 
original matrix; in addition, the reduced matrix equation is 
small and well-conditioned. Mutual coupling between the 
CBFs are accounted for efficiently by using the Secondary 
CBFs in combination with a thresholding scheme to retain 
only the significant ones. The efficiency and accuracy of the 
CBFM is demonstrated by applying it to several microstrip 
example problems, and by comparing the resulting circuit 
parameters with those obtained by using the conventional 
MoM. 

Next, we have presented an efficient algorithm, called the 
CBFM-EMA, for rapid prototyping of microstrip circuits and 
antennas etched on layered media. The CBFM-EMA, which 
can be employed for a preliminary design, replaces the 
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stratified medium with a semi-infinite homogeneous 
grounded slab whose Green’s Functions can be derived 
analytically. 


The method has shown to provide a far better accuracy 
than the Equivalent Circuit approach, and is more general as 
well. The solutions obtained by using the EMA algorithm 
agree well with those generated by using the conventional 
MoM, and the solution process offers considerable time 
advantages over the MoM-only approach. Furthermore, 
almost always the CBFM-EMA results have been found to 
be closer to those obtained from the application of the 
conventional MoM, than are the FEM-derived results. The 
relative time advantage of the EMA makes it a desirable 
choice for rapid prototyping of planar antennas and circuits 
encountered in wireless applications. 


In the second part of paper, the CBFM approach has been 
generalized to the analysis of three dimensional scattering 
problems. The proposed technique bypasses the computation 
of the Secondary CBFs, enabling us to evaluate the reduced 
matrix by using only the Primary CBFs generated by 
illuminating the block with a spectrum of plane waves. In 
particular, we have introduced a versatile and numerically 
efficient computational technique to model the problem of 
scattering from plasmonic nano-rod antennas. The key to 
achieving the numerical efficiency is to utilize macro basis 
functions that take into account the physics of the problem to 
reduce the size of matrix equation to be solved. Next, the 
developed formulation is generalized to characterize the 
reflection and transmission coefficients of a periodic array of 
plasmonic nano-rods illuminated by an obliquely incident 
plane wave. The proposed approach, which uses the concept 
of progressively increasing rings and the Parseval's theorem, 
yields very accurate results and is characterized by a sensible 
time-saving with respect to the conventional MoM approach. 
Finally, we have analyzed large, non-uniform optical antenna 
arrays composed of nano-rods which form optical vortices, 
axicons and lenses by using the CBFM. The CBFM not only 
reduces the number of unknowns significantly without 
scarifying the final accuracy, but also enables us to handle 
large, truncated and non-uniform arrays taking into account 
all the mutual coupling effects in a rigorous manner. 
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